* schneider-troeger-makeplots.do
* This file uses the simulated beta coefficients to make our plots. also does delta method
* -----------------------------------------------------------------------------
cd "~/Dropbox/Benton-Jordan-Philips/final-JOP-replication"
clear
clear matrix
clear mata
set maxvar 80000

* some globals to simplify code below:
global myreshape "reshape long sigma2true_ sigma2_ll95_p sigma2_ll90_p sigma2_ll75_p sigma2_med_p sigma2_ul75_p sigma2_ul90_p sigma2_ul95_p sigma2_ll95_m sigma2_ll90_m sigma2_ll75_m sigma2_med_m sigma2_ul75_m sigma2_ul90_m sigma2_ul95_m sigma2_ll95_r sigma2_ll90_r sigma2_ll75_r sigma2_med_r sigma2_ul75_r sigma2_ul90_r sigma2_ul95_r sigma2_SDul_p sigma2_SDll_p sigma2_SDul_m sigma2_SDll_m sigma2_SDul_r sigma2_SDll_r sigma2_pred_d sigma2_ll95_d sigma2_ul95_d, j(time) i(count)"
global graphhome "Bootstrap-Simulation/fall-2023-RR-submission/figures"


* -----------------
* Create program to re-load sims for each scenario to avoid using too many preserve/restore commands:
cap prog drop reload
prog define reload 
import excel "data/schneider-troeger/st-data.xlsx", sheet("Sheet1")  firstrow clear
* everythign allready sorted
gen t = _n
tsset t
keep diff_dji golf_severity golf_possum golf_negsum pal_severity pal_possum pal_negsum yug_severity yug_possum yug_negsum diff_cac diff_ftse y91 y92 y93 y94 y95 y96 y97 y98 y99 y00 t


*they include a LDV:
gen ldiff_dji = l.diff_dji
compress
set seed 03858420

* merge in all the betas across the 3 approaches:
gen bootcount = _n // to merge on
joinby bootcount using "data/schneider-troeger/schneider-troeger-betas-parametric.dta", unmatched(both) // need unmatched or else it'll cut to T=1000
forv i = 1/30 {
	rename b`i' b_p`i'
}
drop _merge
joinby bootcount using "data/schneider-troeger/schneider-troeger-betas-me.dta", unmatched(both) // need unmatched or else it'll cut to T=1000
forv i = 1/30 {
	rename b`i' b_m`i'
}

drop _merge
joinby bootcount using "data/schneider-troeger/schneider-troeger-betas-resid.dta", unmatched(both) // need unmatched or else it'll cut to T=1000
forv i = 1/30 {
	rename b`i' b_r`i'
}
drop _merge
* b_p = parametric, b_m = maximum entropy, b_r = residual (using enders' starting values)
* We'll create 6 measures of uncertainty: percentiles and +/- SD*1.96 (2 each across the 3 approaches). The line will be the estimated sigma^2 from the OG model for SD method, and median for percentiles


/* 
	Here's the coding of each beta:
	b1 =  diff_dji:golf_sever~y
	b2 =  diff_dji:golf_possum
	b3 =  diff_dji:golf_negsum
	b4 =  diff_dji:pal_severity
	b5 =  diff_dji:pal_possum
	b6 =  diff_dji:pal_negsum
	b7 =  diff_dji:yug_severity
	b8 =  diff_dji:yug_possum
	b9 =  diff_dji:yug_negsum
	b10 = diff_dji:l.diff_dji
	b11 = diff_dji:diff_cac
	b12 = diff_dji:diff_ftse
	b13 = diff_dji:_cons 
	b14 = HET: golf_sever~y
	b15 = HET: pal_sever~y
	b16 = HET: yug_severity
	b17 = HET: y91
	b18 = HET: y92
	b19 = HET: y93
	b20 = HET: y94
	b21 = HET: y95
	b22 = HET: y96
	b23 = HET: y97
	b24 = HET: y98
	b25 = HET: y99
	b26 = HET: y00
	b27 = HET: _cons
	b28 = HET: L.arch
	b29 = HET: L.tarch
	b30 = HET: L.garch
*/

* here's the model
arch diff_dji golf_severity golf_possum golf_negsum pal_severity pal_possum pal_negsum yug_severity yug_possum yug_negsum ldiff_dji diff_cac diff_ftse, arch(1) tarch(1) garch(1) het(golf_severity pal_severity yug_severity y91 y92 y93 y94 y95 y96 y97 y98 y99 y00) arch0(xb)
mat B = e(b)
mat list B
* save true B's for type 2 CI creation
svmat B, names(trueB_b)
forv i = 1/30 {
	su trueB_b`i' in 1
	replace trueB_b`i' = r(mean)
}

* so we care about b14/b30
* assume we're in 1995 (mid-sample). meaning b17/b26 = 0, but b21 = 1j. All other severity variables = 0
end


* -------------------------------------------------------------------------
*  Table 3 (scroll up to GARCH model)
* -------------------------------------------------------------------------
reload // re-load sims for this scenario
* -------------------------------------------------------------------------

* -------------------------------------------------------------------------
*  Figure SM.4
* -------------------------------------------------------------------------
twoway kdensity b_p30 || kdensity b_m30 || kdensity b_r30, xline(.6964717, lcolor(black) lwidth(medthick) lpattern(solid)) xline(.6180526, lcolor(black) lwidth(medthin) lpattern(dash)) xline(.7748908, lcolor(black) lwidth(medthin) lpattern(dash)) legend(order(1 "Parametric" 2 "Maximum Entropy" 3 "Residual") rows(1)) xtitle("Estimated Parameter Value")
*graph export "$graphhome/GARCH_density2.pdf", as(pdf) replace

* -------------------------------------------------------------------------
*  Figure SM.5
* -------------------------------------------------------------------------
twoway kdensity b_p15 || kdensity b_m15 || kdensity b_r15, xline(.3067963, lcolor(black) lwidth(medthick) lpattern(solid)) xline(.2559787, lcolor(black) lwidth(medthin) lpattern(dash)) xline(.3576139, lcolor(black) lwidth(medthin) lpattern(dash)) legend(order(1 "Parametric" 2 "Maximum Entropy" 3 "Residual") rows(1)) xtitle("Estimated Parameter Value")
*graph export "$graphhome/diff_ftse_density2.pdf", as(pdf) replace

* -------------------------------------------------------------------------
*  Figure SM.6
* -------------------------------------------------------------------------
twoway kdensity b_p15 || kdensity b_m15 || kdensity b_r15, xline(.3185028, lcolor(black) lwidth(medthick) lpattern(solid)) xline(.040412, lcolor(black) lwidth(medthin) lpattern(dash)) xline(.5965936, lcolor(black) lwidth(medthin) lpattern(dash)) legend(order(1 "Parametric" 2 "Maximum Entropy" 3 "Residual" ) rows(1)) xtitle("Estimated Parameter Value")
*graph export "$graphhome/pal_severity_density2.pdf", as(pdf) replace


* ----------------------------------------------------------------------------
*	SCENARIO 1: Palestininan-Israeli shock (+1points) 
* ----------------------------------------------------------------------------
* first create OG predictions:
gen sigma2true_0 = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*0 + trueB_b16*0) // golf, pal, and yug severity, respectively
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*0 + trueB_b16*0) + /// golf, pal, and yug severity, respectively 
	trueB_b28*0 + trueB_b29*0 + trueB_b30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
}
* Shock: at t = 31, shock pal by 1 for 1 period:
su sigma2true_30, meanonly
gen sigma2true_31 = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*1 + trueB_b16*0) + /// golf, pal, and yug severity, respectively 
	trueB_b28*0 + trueB_b29*0 + trueB_b30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
* post-shock, at t = 32+, set shock pal back to 0 and keep predicting
forv i = 32/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*0 + trueB_b16*0) + /// golf, pal, and yug severity, respectively 
	trueB_b28*0 + trueB_b29*0 + trueB_b30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
* init (all set to 0, except constant and 95 dummy)
nlcom exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*0 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0) 
mat delta_b = r(table)
* loop, t=1/30
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + _b[ARCH:L.arch]*0 + ///
	_b[ARCH:L.tarch]*0 + ///
	exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*0 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0)
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t= 31, shock pal by 1 for 1 period
nlcom _b[ARCH:L.garch]*`delta_b' + _b[ARCH:L.arch]*0 + ///
	_b[ARCH:L.tarch]*0 + ///
	exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*1 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0)
	mat delta_b = r(table)
	gen sigma2_pred_d31 = delta_b[1,1] // expected value
	gen sigma2_ll95_d31 = delta_b[5,1] // ll
	gen sigma2_ul95_d31 = delta_b[6,1] // ul
* post-shock, at t=32+, pal back to 0 and keep predicting
forv i = 32/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + _b[ARCH:L.arch]*0 + ///
	_b[ARCH:L.tarch]*0 + ///
	exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*0 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0)
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}	
	

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r
	gen sigma2_`b'0 = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*0 + b_`b'16*0) // golf, pal, and yug severity, respectively
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*0 + b_`b'16*0) + /// golf, pal, and yug severity, respectively 
	b_`b'28*0 + b_`b'29*0 + b_`b'30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
	}
	* shock: at t=31, shock pal severity by 1:
	su sigma2_`b'30, meanonly
	gen sigma2_`b'31 = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*1 + b_`b'16*0) + /// golf, pal, and yug severity, respectively 
	b_`b'28*0 + b_`b'29*0 + b_`b'30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2

	* post-shock, at t = 32+, keep tradeval at + 4 and keep predicting
	forv i = 32/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*0 + b_`b'16*0) + /// golf, pal, and yug severity, respectively 
	b_`b'28*0 + b_`b'29*0 + b_`b'30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
	keep in 29/40
	drop time
	gen time = _n - 1
	
* -----------------------------------------------------------------------------
*	Figure SM.17 (a)
* -----------------------------------------------------------------------------
	* Parametric percentile
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/ST_palshock-parametric-percentile.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.17 (b)
* -----------------------------------------------------------------------------
	* ME percentile
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/ST_palshock-me-percentile.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.17 (c)
* -----------------------------------------------------------------------------
	* resid percentile
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are residual bootstrap percentiles, residual bootstrap approach")

* -----------------------------------------------------------------------------
*	Figure SM.18 (a)
* -----------------------------------------------------------------------------	
	* now parametric SD
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/ST_palshock-parametric-sd.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.18 (b)
* -----------------------------------------------------------------------------
	* ME SD
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/ST_palshock-me-sd.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure SM.18 (c)
* -----------------------------------------------------------------------------	
	* resid SD
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are +/- 1.96*SD residual bootstrap")
	*graph export "$graphhome/ST_palshock-resid-sd.pdf", as(pdf) replace

	* NOT SHOWN
	* delta (for delta we'll show the underlying delta expectation)
	*twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs and expected values from delta method")
	*graph export "$graphhome/ST_palshock-delta-sd.pdf", as(pdf) replace

* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}


* -----------------------------------------------------------------------------
*	Figure 4 (a)
* -----------------------------------------------------------------------------
	* parametric
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/ST_palshock-parametric-percentile-rescale.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure 4 (b)
* -----------------------------------------------------------------------------
	* ME percentile
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/ST_palshock-me-percentile-rescale.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure 4 (c)
* -----------------------------------------------------------------------------
	* resid percentile
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs are residual bootstrap percentiles, residual bootstrap approach")
	*graph export "$graphhome/ST_palshock-resid-percentile-rescale.pdf", as(pdf) replace
	
	* NOT SHOWN
	* delta:
	*twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method")
	*graph export "$graphhome/ST_palshock-delta-percentile-rescale.pdf", as(pdf) replace
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad


* -----------------------------------------------------------------------------
*	Figure SM.19 (a)
* -----------------------------------------------------------------------------
* parametric
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/ST_palshock-parametric-effectsplot.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.19 (b)
* -----------------------------------------------------------------------------
* ME
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/ST_palshock-me-effectsplot.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.19 (c)
* -----------------------------------------------------------------------------
* resid
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/ST_palshock-resid-effectsplot.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure SM.19 (d)
* -----------------------------------------------------------------------------
* Delta
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method")
	*graph export "$graphhome/ST_palshock-delta-effectsplot.pdf", as(pdf) replace
* ----------------------------------------------------------------------------





reload // re-load sims for this scenario




* ----------------------------------------------------------------------------
*	Scenario 2: Palestininan-Israeli shock (+1points), ongoing for 3 days
* ----------------------------------------------------------------------------
* first create OG predictions:
gen sigma2true_0 = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*0 + trueB_b16*0) // golf, pal, and yug severity, respectively
forv i = 1/30 {
	loc pasti = `i' - 1
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*0 + trueB_b16*0) + /// golf, pal, and yug severity, respectively 
	trueB_b28*0 + trueB_b29*0 + trueB_b30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
}
* shock: at t=31, 32 and 33, shock pal severity by 1:
forv i = 31/33 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*1 + trueB_b16*0) + /// golf, pal, and yug severity, respectively 
	trueB_b28*0 + trueB_b29*0 + trueB_b30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
}
* post-shock, at t = 34+, set shock pal back to 0 and keep predicting
forv i = 34/50 {
	su sigma2true_`pasti', meanonly
	gen sigma2true_`i' = exp(trueB_b27 + trueB_b17*0 + trueB_b18*0 + trueB_b19*0 + trueB_b20*0 + trueB_b21*1 + trueB_b22*0 + trueB_b23*0 + trueB_b24*0 + trueB_b25*0 + trueB_b26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	trueB_b14*0 + trueB_b15*0 + trueB_b16*0) + /// golf, pal, and yug severity, respectively 
	trueB_b28*0 + trueB_b29*0 + trueB_b30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
}

* now for delta (saving as sigma2_ll95_d, sigma2_ul95_d, sigma2_pred_d)
* init (all set to 0, except constant and 95 dummy)
nlcom exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*0 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0) 
mat delta_b = r(table)
* loop, t=1/30
forv i = 1/30 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + _b[ARCH:L.arch]*0 + ///
	_b[ARCH:L.tarch]*0 + ///
	exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*0 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0)
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* shock: at t= 31/33, shock pal by 1 for 1 period
forv i = 31/33 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + _b[ARCH:L.arch]*0 + ///
	_b[ARCH:L.tarch]*0 + ///
	exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*1 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0)
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}
* post-shock, at t=34+, pal back to 0 and keep predicting
forv i = 34/50 {
	loc delta_b = delta_b[1,1]
	nlcom _b[ARCH:L.garch]*`delta_b' + _b[ARCH:L.arch]*0 + ///
	_b[ARCH:L.tarch]*0 + ///
	exp(_b[HET:_cons] + _b[HET:golf_severity]*0 + _b[HET:pal_severity]*0 + ///
	_b[HET:yug_severity]*0 + _b[HET:y91]*0 + _b[HET:y92]*0 + _b[HET:y93]*0 + ///
	_b[HET:y94]*0 + _b[HET:y95]*1 + _b[HET:y96]*0 + _b[HET:y97]*0 + ///
	_b[HET:y98]*0 + _b[HET:y99]*0 + _b[HET:y00]*0)
	mat delta_b = r(table)
	gen sigma2_pred_d`i' = delta_b[1,1] // expected value
	gen sigma2_ll95_d`i' = delta_b[5,1] // ll
	gen sigma2_ul95_d`i' = delta_b[6,1] // ul
}

* now do the same for the bootstrap procedures:
foreach b in p m r { // for b_p b_m b_r
	gen sigma2_`b'0 = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*0 + b_`b'16*0) // golf, pal, and yug severity, respectively
	forv i = 1/30 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*0 + b_`b'16*0) + /// golf, pal, and yug severity, respectively 
	b_`b'28*0 + b_`b'29*0 + b_`b'30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
	}
	* shock: at t=31, 32 and 33, shock pal severity by 1:
	forv i = 31/33 {
		loc pasti = `i' - 1
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*1 + b_`b'16*0) + /// golf, pal, and yug severity, respectively 
	b_`b'28*0 + b_`b'29*0 + b_`b'30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
	}
	
	* post-shock, at t = 34+, keep tradeval at + 4 and keep predicting
	forv i = 34/50 {
		su sigma2_`b'`pasti', meanonly
		gen sigma2_`b'`i' = exp(b_`b'27 + b_`b'17*0 + b_`b'18*0 + b_`b'19*0 + b_`b'20*0 + b_`b'21*1 + b_`b'22*0 + b_`b'23*0 + b_`b'24*0 + b_`b'25*0 + b_`b'26*0 + /// constant (b27) and year dummies set to zero, except 1995 (b21)
	b_`b'14*0 + b_`b'15*0 + b_`b'16*0) + /// golf, pal, and yug severity, respectively 
	b_`b'28*0 + b_`b'29*0 + b_`b'30*r(mean) // arch terms = 0, garch term (b30) = l.sigma2
	}
}
* create CIs:
foreach b in p m r {
	forv i = 1/50 {
		* Type 1: %iles:
		_pctile sigma2_`b'`i', p(2.5, 5, 12.5, 50, 87.5, 95, 97.5) 
		gen sigma2_ll95_`b'`i' = r(r1)
		gen sigma2_ll90_`b'`i' = r(r2)
		gen sigma2_ll75_`b'`i' = r(r3)
		gen sigma2_med_`b'`i' = r(r4)
		gen sigma2_ul75_`b'`i' = r(r5)
		gen sigma2_ul90_`b'`i' = r(r6)
		gen sigma2_ul95_`b'`i' = r(r7)
		
		* Type 2: OG prediction +/- 1.96*SD
		su sigma2_`b'`i'
		gen sigma2_SDul_`b'`i' =  sigma2true_`i' + 1.96*r(sd)
		gen sigma2_SDll_`b'`i' =  sigma2true_`i' - 1.96*r(sd)
	}
}

* now collapse, reshape, drop burn-ins and plot:
keep sigma2*
drop sigma2true_0 // want time to start at 1
collapse sigma2*
gen count = _n
$myreshape
	keep in 29/40
	drop time
	gen time = _n - 1
	

* -----------------------------------------------------------------------------
*	Figure SM.20 (a)
* -----------------------------------------------------------------------------
	* Parametric percentile
	twoway rarea sigma2_ll95_p sigma2_ul95_p time, color(eltblue) || rarea sigma2_ll90_p sigma2_ul90_p time, color(ebblue) || rarea sigma2_ll75_p sigma2_ul75_p time, color(edkblue) || line sigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-parametric-percentile.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.20 (b)
* -----------------------------------------------------------------------------
	* ME percentile
	twoway rarea sigma2_ll95_m sigma2_ul95_m time, color(eltblue) || rarea sigma2_ll90_m sigma2_ul90_m time, color(ebblue) || rarea sigma2_ll75_m sigma2_ul75_m time, color(edkblue) || line sigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-me-percentile.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.20 (c)
* -----------------------------------------------------------------------------
	* resid percentile
	twoway rarea sigma2_ll95_r sigma2_ul95_r time, color(eltblue) || rarea sigma2_ll90_r sigma2_ul90_r time, color(ebblue) || rarea sigma2_ll75_r sigma2_ul75_r time, color(edkblue) || line sigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-resid-percentile.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure SM.21 (a)
* -----------------------------------------------------------------------------
	* parametric SD
	twoway rarea sigma2_SDll_p sigma2_SDul_p time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are +/- 1.96*SD parametric bootstrap")
	*graph export "$graphhome/ST_pal3dayshock-parametric-sd.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.21 (b)
* -----------------------------------------------------------------------------
	* ME SD
	twoway rarea sigma2_SDll_m sigma2_SDul_m time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are +/- 1.96*SD maximum entropy bootstrap")
	*graph export "$graphhome/ST_pal3dayshock-me-sd.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.21 (c)
* -----------------------------------------------------------------------------
	* resid SD
	twoway rarea sigma2_SDll_r sigma2_SDul_r time, color(eltblue) || line sigma2true_ time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs are +/- 1.96*SD residual bootstrap")
	graph save g4.gph, replace
	*graph export "$graphhome/ST_pal3dayshock-resid-sd.pdf", as(pdf) replace
	
	* NOT SHOWN
	* delta (for delta we'll show the underlying delta expectation)
	*twoway rarea sigma2_ll95_d sigma2_ul95_d time, color(eltblue) || line sigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance") note("CIs and expected values from delta method")
	

* now express as % of pre-shock value:
foreach b in p m r {
	su sigma2_med_`b' if time == 1
	foreach var of varlist sigma2_med_`b' sigma2_ll95_`b' sigma2_ul95_`b' sigma2_ll90_`b' sigma2_ul90_`b' sigma2_ll75_`b' sigma2_ul75_`b' {
		gen alt`var' = (`var'*100)/r(mean)
	}
}
* do the same for delta:
su sigma2_pred_d if time == 1
foreach var of varlist sigma2_pred_d sigma2_ll95_d sigma2_ul95_d {
	gen alt`var' = (`var'*100)/r(mean)
}

* -----------------------------------------------------------------------------
*	Figure 5 (a)
* -----------------------------------------------------------------------------
* Parametric percentile
	twoway rarea altsigma2_ll95_p altsigma2_ul95_p time, color(eltblue) || rarea altsigma2_ll90_p altsigma2_ul90_p time, color(ebblue) || rarea altsigma2_ll75_p altsigma2_ul75_p time, color(edkblue) || line altsigma2_med_p time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-parametric-percentile-rescale.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure 5 (b)
* -----------------------------------------------------------------------------	
	* ME percentile
	twoway rarea altsigma2_ll95_m altsigma2_ul95_m time, color(eltblue) || rarea altsigma2_ll90_m altsigma2_ul90_m time, color(ebblue) || rarea altsigma2_ll75_m altsigma2_ul75_m time, color(edkblue) || line altsigma2_med_m time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-me-percentile-rescale.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure 5 (c)
* -----------------------------------------------------------------------------
	* resid percentile
	twoway rarea altsigma2_ll95_r altsigma2_ul95_r time, color(eltblue) || rarea altsigma2_ll90_r altsigma2_ul90_r time, color(ebblue) || rarea altsigma2_ll75_r altsigma2_ul75_r time, color(edkblue) || line altsigma2_med_r time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-resid-percentile-rescale.pdf", as(pdf) replace
	
	* NOT SHOWN
	* delta:
	* twoway rarea altsigma2_ll95_d altsigma2_ul95_d time, color(eltblue) || line altsigma2_pred_d time, lcolor(black) legend(off) ytitle("Expected Conditional Error Variance" "(as % of Pre-Shock Variance)") note("CIs and expected values from delta method")
	
* Now create effectsplot-style interpretation:
foreach b in p m r { // for bootstraps
	foreach z in med ll95 ll90 ll75 ul95 ul90 ul75 {
		su sigma2_med_`b' if time == 1 // baseline time (i.e., just before shock)
		gen sigma2_`z'_`b'_sr = sigma2_`z'_`b' - r(mean) 
		gen sigma2_`z'_`b'_lr = sigma2_`z'_`b'_sr
		replace sigma2_`z'_`b'_sr = . if time != 2 // this is the SR effect
		replace sigma2_`z'_`b'_lr = . if time != 11 // this is the LR effect
	}
}
foreach var of varlist sigma2_ll95_d sigma2_ul95_d sigma2_pred_d { // for delta
	su sigma2_pred_d if time == 1 // baseline time (just before shock)
	gen `var'_sr = `var' - r(mean)
	gen `var'_lr = `var'_sr
	replace `var'_sr = . if time != 2 // this is the SR effect
	replace `var'_lr = . if time != 11 // this is the LR effect
}
* now plot:
gen sortplot = .
replace sortplot = 1 in 3 // SR
replace sortplot = 1.5 in 12 // LR
replace sortplot = .5 in 1 // pad
replace sortplot = 2 in 2 // pad


* -----------------------------------------------------------------------------
*	Figure SM.22 (a)
* -----------------------------------------------------------------------------
* parametric
twoway rspike sigma2_ll95_p_sr sigma2_ul95_p_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_p_lr sigma2_ul95_p_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_p_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_p_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are parametric bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-parametric-effectsplot.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure SM.22 (b)
* -----------------------------------------------------------------------------
* ME
twoway rspike sigma2_ll95_m_sr sigma2_ul95_m_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_m_lr sigma2_ul95_m_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_m_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_m_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are maximum entropy bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-me-effectsplot.pdf", as(pdf) replace

* -----------------------------------------------------------------------------
*	Figure SM.22 (c)
* -----------------------------------------------------------------------------
* resid
twoway rspike sigma2_ll95_r_sr sigma2_ul95_r_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_r_lr sigma2_ul95_r_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_med_r_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_med_r_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs are residual bootstrap percentiles")
	*graph export "$graphhome/ST_pal3dayshock-resid-effectsplot.pdf", as(pdf) replace
	
* -----------------------------------------------------------------------------
*	Figure SM.22 (d)
* -----------------------------------------------------------------------------
* Delta
twoway rspike sigma2_ll95_d_sr sigma2_ul95_d_sr sortplot, horizontal ///
	lcolor(black) || rspike sigma2_ll95_d_lr sigma2_ul95_d_lr sortplot, ///
	horizontal lcolor(black) || scatter sortplot sigma2_pred_d_sr, ///
	msymbol(T) mcolor(black) msize(medlarge) || scatter sortplot ///
	sigma2_pred_d_lr, msymbol(S) mcolor(black) msize(medlarge) ///
	xline(0, lcolor(black) lstyle(solid)) yscale(axis(1) noline) ///
	xlabel(, grid glcolor(gs15)) plotregion(style(none)) ///
	ylabel(1 "Short-Run" 1.5 "Long-Run") legend(off) ///
	xtitle("Expected Change from Pre-Shock Variance") ytitle("") ///
	note("CIs and expected values from delta method")
	*graph export "$graphhome/ST_pal3dayshock-delta-effectsplot.pdf", as(pdf) replace

* ----------------------------------------------------------------------------
